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Preface 

This  report  is  the  written  version  of 
a  talk  presented  by  William  Grossmann 
at  the  Varenna  Workshop  on  High  Beta 
Plasmas  at  Varenna,  Italy,  September 
1977.   It  presents  some  of  the  high- 
lights of  the  extensive  theoretical 
and  numerical  program  on  plasma  equi- 
librium, transport,  and  adiabatic 
processes,  including  an  outline  of 
some  of  the  simulation  codes  which 
have  been  developed  to  implement 
this  body  of  work. 
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Abstract 

A  survey  is  given  of  the  development  of 
classical  diffusion  theory  which  arose  from 
the  observation  of  Grad  and  Hogan  that  the 
Pf irsch-Schliiter  and  Neoclassical  theories 
are  very  special  and  frequently  inapplic- 
able because  they  require  that  plasma  mass 
flow  be  treated  as  transport  rather  than  as 
a  state  variable  of  the  plasma.   The  sub- 
sequent theory,  efficient  numerical 
algorithms,  and  results  of  various  operating 
codes  are  described. 


The  theory  of  classical  diffusion  in  toroidal  (or 
periodic)  systems  originated  with  Kruskal  and  Kulsrud  in 
1958  [1].   They  took  a  form  of  Ohm's  law  equivalent  to 

(1)  E  +  u  X  B  =  nj 

and  assumed  that  the  system  was  fully  stationary  in  time 
[pdu/dt  dropped  from  momentum  conservation ,  3B/9t  from  flux 
conservation,  and  8p/3t  from  mass  conservation  -  except  that 
an  artificial  mass  source  was  inserted  in  order  to  avoid  the 
consequent  restriction  to  only  trivial  solutions] .   Pfirsch 
and  Schliiter  [2]  made  this  theory  more  useful  by  evaluating 
explicit  formulas  for  the  diffusion  rate  in  the  special  case 
of  large  aspect  and  approximately  circular  flux  surfaces. 
Some  time  later  these  formulas  were  generalized  by  Maschke  [3] 
to  include  c|>  E*dx  7^  0  (consistent  with  curl  E  =  0  and 
9B/3t  =  0  in  the  plasma  itself)  in  order  to  represent  the 
Tokamak  configuration. 

There  are  two  outstanding  mysteries  which  arise  in  this 
"classical"  theory;  these  were  pointed  out  (and  resolved)  in 
1970  by  Grad  and  Hogan  [4]  and  Grad  [5].   A  velocity  field 
u(x,y,z)  is  set  up  by  the  resistivity  (in  what  would  be  a 
static  equilibrium,  ideally)  but  only  the  net  flow  through 
a  flux  surface 


(2)       U  =  (t  u'ds  =  <u'Vp>/p'  (V)  =  <u*Vv>   ,   V  =  Volume 
is  of  major  interest.   Surprise  number  one  is  that  the  net 


flow  can  be  evaluated  directly  without  having  to  first 
compute  the  entire  flow  field  u(x,y,z).   Multiplying  (1) 
by  J  and  averaging  over  a  flux  surface, 

2 

(3)  <E'J>  -  <U'J  X  B>  =  <nJ  >  • 

By  assumption 

(4)  Vp  =  J  X  B 

(5)  curl  E  =  0 

from  which  the  first  term  in  (3)  vanishes  and  we  obtain 

(6)  U  =  <nJ^>/p' (V) 

which  is,  in  principle,  computable  in  terms  of  a  given  static 

equilibrium.   If  J'B  =  0,  then  J  =  i ^P I  /B  ,  and  the 

2 
classical  1/B   formula  results 


J) 


<^^>- 


The  extra  contribution  from  J  .  with  approximately  circular 
cross-section  gives  the  familiar  Pf irsch-Schliiter  factor. 

The  second,  more  disturbing  consequence  of  the  Kruskal , 
Kulsrud,  Pfirsch,  Schliiter  (K  PS)  theory  is  the  fact  that  this 
model  severely  restricts  the  relevant  equilibria.   In  a 
general  toroidal  configuration  (ignoring  the  very  real  questions 
of  a  nonexistence  of  equilibria  through  topological  instability, 
Grad  1967  [6]),  one  should  be  able  to  specify  two  essentially 


arbitrary  profiles  (cf.  Grad  and  Rubin  [7]  or  Kruskal  and 

2 
Kulsrud  [1] ) .   The  K  PS  diffusion  theory  and  formulas  apply 

only  to  a  special  subclass  of  equilibria  for  which 

(8)  <nJ'B>   =   0 

on  each  surface  (an  immediate  consequence  of  (1) ,  after 
multiplication  by  B  and  averaging) .   For  the  Tokamak, 
o  E'dx  =  c  7^  0,  the  restriction 

(9)  C4;'  (V)  =  <nJ-B> 

also  allows  only  one  profile  to  be  freely  specified. 
The  two  questions, 

i)  Why  is  the  profile  restricted,  and  what  is  the 
diffusion  rate  in  general,  without  this  restric- 
tion? 

ii)  What  is  the  source  of  the  miracle,  that  the  net 
flow  can  be  calculated  without  first  solving  a 
boundary  value  problem  for  u(x,y,z),  and  does 
it  persist  in  the  general  case? 

are  resolved  by  dropping  curl  E  =  0  and  re-introducing  3B/8t. 
The  extension  to  curl  E  7^  0  turns  out  to  allow  general 
equilibriiam  profiles,  gives  entirely  different  "paleo- 
classical"  diffusion  results,  and  requires  solution  of  a 
global  boundary  value  problem  before  obtaining  the  desired 
surface  averages.   It  also  introduces  coupling  between  plasma 
diffusion  and  another  classical  phenomenon,  the  skin  effect. 


Question  (ii)  above  is  answered  in  greater  generality  in 
[5] .   Briefly,  in  an  irreversible  thermodynamic  formulation 

(10)  V.  =  A. .  u. 

1     ID  3 

in  which  the  fluxes  v. (x) ,  forces  u.(x)  and  transport  co- 
efficients A. . (x)  depend  on  a  parameter  x.   The  net  flux 
/  v.dx  can  be  a  linear  function  of  the  forces  with  a  universal 
matrix  which  is  independent  of  the  values  of  u-  only  when 
either  the  u.  or  v.  are  independent  of  x  (in  which  case  the 
simplification  is  evident) .   In  the  plasma  diffusion  example, 
the  facts  that  the  "force"  p  [more  accurately,  p'(V)]  is 
constant  on  the  flux  surface,  and  the  force  (actually  flux!) 
E*dx  is  path  independent,  are  what  allow  surface  averaged 
transport  coefficients.   Similarly,  with  heat  flow  the 
assumption  that  T  is  constant  on  a  flux  surface  must  be  made. 
Without  these  assiamptions  (e.g.  with  anisotropic  stress 
tensor,  or  T  7^  T(ip),  or  curl  E  7^  0)  there  are  no  global  trans- 
port coefficients  relating  surface  averaged  fluxes;  the  entire 
time-dependent,  two  or  three  dimensional  flow  problem  must  be 
solved,  after  which  averages  can  be  computed;  (however,  by 
various  ingenious  stratagems  and  several  years  of  development, 
the  total  amount  of  effort  required  to  solve  the  more  general 
transport  problems  is  not  much  more  than  for  the  restricted 
K  PS  special  case) . 

The  mathematical  and  numerical  problem  after  reinsertion 
of  9B/9t  is  very  nonstandard.   The  reason  is  found  by  comparing 
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this  Grad-Hogan  (GH)  model  with  the  full  MHD  system  including 
Pdu/dt  as  well  as  3B/8t.   Instead  of  marching  u  in  time,  the 
GH  system  states  that  u  must  be  determined  by  the  constraint 
that  B  and  p  continue  to  satisfy  J  x  B  =  Vp  as  they  vary  in 
time.   The  canonical  procedure  for  obtaining  the  equation 
governing  u  is  to  differentiate  the  constraint  (i.e.  the 
pressure  balance)  with  respect  to  t,  replacing  3p/8t,  3B/8t, 
3j/9t  with  expressions  involving  u.   The  resulting  monstrous 
equation  for  u  is  nonstandard,  and  is  not  a  differential 
equation.   In  2-D,  axial  symmetry,  and  helical  symmetry  it  is 
a  linear  GDE  (originally  QDE) ,  for  which  substantial  theory 
and  effective  numerical  algorithms  have  been  developed  since 
1970  (theory  [8-10],  numerical  [11-16]).   There  is  a  com- 
putational tradeoff  between  taking  a  conventional  system  of 
differential  equations,  including  pdu/dt,  with  its  very  wide 
span  of  physical  time  scales,  or  taking  the  less  conventional 
GDE  on  a  much  more  manageable  transport  time  scale.   Because 
of  the  substantial  analytical  and  numerical  developments, 
Refs.  [8-18] ,  the  correct  decision  would  seem  to  be  in  favor 
of  GH  rather  than  full  MHD  with  pdu/dt,  except  for  problems 
of  plasma  initiation  (usually  on  microsecond  time  scales) . 

The  key  to  developing  accurate,  fast,  efficient  plasma 
transport  codes  turned  out  to  depend  on  techniques  developed 
for  the  adiabatic  problem.   The  two  problems  are  mathematically 
similar  in  that  each  follows  a  family  of  static  equilibria. 
The  very  successful  adiabatic  algorithm  depends  on  alternating 


between  1-D  profile  determination  (from  an  ODE  with  tentative 
inductance  coefficients)  and  2-D  contour  determination  [from 
an  elliptic  PDE  with  tentative  current  profile  J (V) ] .   The 
contour  shapes  determine  the  inductance  profiles,  and  the 
average  pressure  balance  (containing  inductance  and  given 
adiabatic  profiles)  determines  the  correct  profile  [9,11]. 
Among  the  important  numerical  features  are  the  use  of 
Aif;  =  J(V)  instead  of  Aijj  =  J  Uj)    to  improve  numerical  stability, 
and  application  of  analytical  knowledge  of  the  singular 
behavior  near  the  plasma  center  (V  =  0) ,  plasma  edge,  and 
separatrix  (different  for  adiabatic  and  resistive,  symmetric 
and  asymmetric  cases,  etc.)  which  allows  much  more  efficient 
computation. 

A  prototype  adiabatic  GDE  is 

(11)  Aip      =    -    ip"  . 

On  the  left  is  an  elliptic  operator  on  ilj(x,y).   Writing  V(4;) 
for  the  volume  (area  in  2-D)  within  a  tJ;-contour,  and  ^  {V) 
for  the  inverse  function,  the  right  side  is  an  ordinary  differ- 
ential operator  on  ^{V).      Conventional  methods  of  solution 
which  work  for  the  elliptic  equation,  A\p   =  f  iii)  ,    or  even 
A^j  =  f(v)  (e.g.  iteration  on  the  right  side)  do  not  work  for 
the  GDE.   The  iteration  between  1-D  and  2-D,  i.e.,  between 

(11)  and  its  average 

(12)  (Kij;')'  =    -    i>"    ,    K(V)  =  <|VV|^> 
is  the  only  known  practical  numerical  method. 


There  are  essentially  three  procedures  for  solving  non- 
inertial  (P  5^=  0)  transport  problems.   The  first  (1970  [8]) 
makes  use  of  the  Green's  function  for  the  elliptic  operator 
which  is  the  linearization  (small  variation  in  time  At)  of 
the  equilibrium.   This  is  very  impractical  numerically  (since 
the  Green's  function  has  to  be  recomputed  each  time  step), 
but  is  very  appropriate  for  analytical  work  and  has  led  to 
an  existence  proof  for  the  GDE  involved  in  this  problem  [10] 
(also  for  the  adiabatic  problem) . 

A  second  procedure  is  to  solve  the  linear  GDE  for  the 
velocity  u'Vijj  at  each  time  step  (solved  by  iterating  between 
profile  and  geometry) ,  inserting  the  velocity  to  then  march 
ii ,      p,  etc.   This  is  numerically  impractical  because  the 
velocity  profile  exhibits  boundary  layers  and  rapid  variation 
and  would  require  small  time  steps.   It  is  occasionally  useful 
analytically,  and  was  the  source  of  the  original  demonstration 
by  Grad  and  Hogan  that  in  some  very  special  cases  (principally 
low  6,  large  aspect,  small  poloidal  field)  there  is  a  rapid 
(skin  effect)  time  scale  which  is  nonclassical  during  which 

which  the  two  independent  equilibrium  profiles  relax  until 

2 
they  satisfy  the  K  PS  profile  restriction,  and  during  which 

the  diffusion  velocity  relaxes  to  the  "classical"  value. 

(A  somewhat  similar  theory  was  repeated  later  [17] ,  and  has 

been  carried  out  with  more  general  transport  [18]). 

The  third  procedure  for  GH  diffusion,  the  one  which  has 

proved  to  be  most  effective  -  so  much  so  that  a  very  nonstandard 
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2-D  diffusion  problem  is  reduced  to  one  with  computation 
time  only  slightly  more  than  for  a  standard  1-D  diffusion 
problem  -  is  to  introduce  independent  and  dependent  variables 
so  as  to  eliminate  the  convection  velocity;  (it  can  be  computed 
afterwards,  if  desired) .   The  time  steps  can  be  taken  very 
large  because  (in  the  proper  variables)  the  geometric  coef- 
ficients vary  very  slowly,  even  through  a  change  in  topology. 
In  a  moderately  optimized  resistive  calculation  in  which  a 
Belt  Pinch  is  squeezed,  islates,  and  develops  into  a  fully- 
formed  Doublet,  the  computation  time  on  a  CDC  6600  was  one 
minute. 

The  crucial  point  is  to  take  the  proper  independent  and 
dependent  variables  in  the  properly  averaged  transport  and 
pressure  balance  equations.   For  example,  ii ,    X/  ^  =  /  pdV 
are  possible  independent  variables;  V  is  not  (except  at  low  3, 
large  aspect,  small  poloidal  field) .   As  dependent  variables, 
^   and  X  are  possible;  ijj'(V)  is  not,  but  ii'/p    is  satisfactory; 
p  is  not  a  good  variable  but  p/p^  or  p/(tj;')   are;  with 
reversed  field,  ijj  or  x  is  a  poor  independent  variable.   More 
precisely,  normalized  independent  variables,  e.g.  V/V^^   or 

\b/\l)    ^  are  essential  because  the  inductance  coefficients 

^'  plasma 

(and  consequently  the  current)  have  moving  singularities, 
following  V    ,  and  a  current  skin  layer  may  follow  V  lagma* 
Before  writing  down  the  mathematical  formulation,  there 
are  several  qualitative  features  to  mention.   The  basic 
concept  is  that  of  a  slowly  varying  family  of  equilibria  (this 
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is  common  to  adiabatic  and  dissipative  models  and  is  the 
cause  of  the  mathematical  and  numerical  similarity).   In  Ohm's 
law,  (1)  J  and  B  can  be  taken  as  given  (they  define  the 
instantaneous  equilibrium) ;  E  and  u  are  to  be  found.   The 
question  is  how  the  emf  nJ  splits  into  its  two  components  E 
and  u  X  B.   If,  for  the  sake  of  argument,  the  split  goes 
entirely  one  way,  E  =  nJ/  the  result  is  a  classical  (but  not 
"classical")  diffusion  equation  for  the  skin  effect. 

(13)  If  "  "  ^^^-^    E  =  -  curl  (n/y^  curl  B) 

with  diffusion  coefficients  D  =  '^/^^q'   "*""  ^^^  opposite  limit, 
u  X  B  =  nJ/ 

u  =  -  nJ  X  B/B^ 
(14) 

=  -  nVp/B^ 

which  is  the  "classical"  diffusion  rate  (assuming,  for  simplicity, 
that  u  is  perpendicular  to  B) .   To  demonstrate  a  diffusion 
equation,  this  velocity  is  inserted  into  mass  conservation  (for 
simplicity,  isothermal) 

(15)  1^  -   div[(np/B^)Vp] 

This  is  a  nonlinear  diffusion  equation  with  diffusion  coefficient 

D,  =  ^  3d  .   The  actual  split  of  nJ  between  E  and  u  x  B  varies 
1    2    o 

from  point  to  point  and  involves  solution  of  a  boundary  value 
problem  (usually  a  GDE) . 
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In   a    low    3,    2-D   model     [8]    the   basic    equation    for   dif- 
fusion  of    flux    is 

(16)  1^      +    u'Vifj   =    nAij; 

where,  from  the  scaling  we  know  that  <u*V4j>  =  0  (higher  order). 
If  we  were  to  ignore  convection, 

(17)  II  =  nA4. 

we  would  be  able  to  solve  for  tjj(x,y,t)  from  this  equation  alone, 
This  solution  would  be  incompatible  with  the  equilibrium 
constraint,  Ai|;  =  f(4^)  {^   contours  and  Ai|;  contours  coincide); 
the  term  u-V^j  is  needed  to  maintain  this  constraint.   The 
average  of  (16)  is 

(18)  ip^      =    n(Ki(j')  • 

2 

where  ij>  (V,t)  is  the  derivative  with  V  fixed  and  K(V)  =  <|VV|  > 

is  a  geometrical  coefficient.   As  diffusion  progresses,  the 
change  in  shape  of  the  iJ;-contours  is  visible  in  (16)  through 
u'ViJ;  and  in  (18)  through  the  changes  in  the  inductance  K(V,t). 
The  second  version,  in  terms  of  K,  is  much  more  intuitive  than 
u'Vijj,  and  it  is  also  more  direct  numerically.   For  example, 
for  a  contour  of  given  circumference,  K  is  a  minimum  for  a 
circle.   It  increases  with  corrugations  or  elongation  of  the 
contour . 

Figure  1  shows  the  history  of  a  Belt  Pinch  squeezed  by 
a  rising  coil  current  at  the  waist  [subject  to  exactly 
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eq.  (16) J.   The  boundary  condition  at  the  edge  is  chosen  to 
eliminate  the  skin  effect.   However,  the  effect  of  corrugation 
of  the  outer  surfaces  is  seen  as  a  rise  in  current  density 
(it  would  be  negative,  as  a  skin  current) .   As  the  point  of 
islation  is  approached,  the  central  contours  elongate, 
increasing  K,  and  causing  a  current  peak.   After  islation,  a 
simple  calculation  shows  that  K  has  a  logarithmic  singularity 
at  the  separatrix.   A  simple  explicit  solution  of  (18)  for  a 
uniformly  moving  K(V-ct)  with  a  logarithmic  singularity  closely 
imitates  the  current  signature  in  Fig.  1.   This  moving  current 
peak  is  seen  in  Doublet  experiments. 

Perhaps  the  most  important  si^ngle  phenomenon  to  be  veri- 
fied by  this  numerical  experiment  is  the  fact  that  current 
penetration  of  a  deformable  medium  (e.g.  plasma)  is  not 
governed  by  the  skin  effect  as  in  a  solid  conductor,  but  the 
current  can  penetrate  instantaneously  (anomalously!)  when 
there  is  a  change  of  shape.   As  a  result  of  a  cyclic  external 
coil  variation,  one  can  even  arrange  to  create  an  AC  (damped 
sinusoidal)  skin  current  at  the  center  of  the  plasma  rather 

than  at  the  edge. 

2 

The  significance  of  the  difference  between  the  K  PS  and 

the  GH  theories  is  sometimes  confused  by  taking  the  techniques 
of  the  latter  (alternating  geometry  and  1-D  transport)  and 
applying  them  to  an  inconsistent  model  which  includes  3B/8t 
and  at  the  same  time  evaluates  mass  flow  by  transport  formulas 
(e.g.  [19]).   It  is  only  in  the  absence  of  8B/8t  (and  with 
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p  =  p(4^),  T  ==  T(i>),    etc.)  that  mass  flow  can  be  described  in 
terms  of  surface  to  surface  transport.   The  overdetermination 
introduced  by  including  8B/3t  and  a  formula  for  mass  trans- 
port is  removed  by  including  only  part  of  Ohm's  law.   There 
is  no  simplification  associated  with  this  ad  hoc  procedure, 
since  the  correct  treatment  of  velocity  (mass  flow)  by  the 

GH  algorithms  is  at  least  as  simple  and  is  more  accurate. 

2 

The  large  discrepance  between  K  PS  "classical"  mass  transport 

and  the  correct  classical  GH  mass  transport  in  the  presence 
of  time  varying  fields  has  been  calculated  [12]  .   The  same 
remarks,  of  course,  apply  to  Neoclassical  transport  which  does 
not  yet  include  transient  fields  correctly  (except  in  some 
special  calculations  such  as  the  Ware  pinch) . 

There  is  also  occasional  misunderstanding  that  the 
use  of  averaged  transport  is  necessarily  a  rough  approximation 
to  the  correct  procedure;  under  appropriate  circumstances  (as 
here  described) ,  it  is  exact  with  no  loss  of  information  as 
demonstrated  by  the  theory  of  the  GDE . 

One  other  qualitative  point  is  the  distinction  between 
islation  and  bifurcation.   The  necessity  of  island  structure  in 
plasma  equilibria  was  pointed  out  in  [6,20,21].   Mathematically, 
bifurcation  (multiple  solutions)  can  be  expected  when  a  point 
eigenvalue  (of  the  dynamical  operator)  crosses  the  origin; 
islation  occurs  (with  unique  continuation  of  the  solution)  when 
a  continuum  touches  the  origin  [20,9].   Each  phenomenon 
requires  numerical  care,  but  of  an  entirely  different  type. 
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The  fact  that  a  change  of  topology  is  compatible  with  ideal 
MHD  [9,22]  increases  the  importance  of  the  adiabatic  model 
(which,  for  example,  gave  the  first  demonstration  of  satura- 
tion of  the  tearing  mode  [9]  and  the  first  time-dependent 
simulation  of  a  topology  change  [12]). 

A  detailed  summary  of  the  equations  and  algorithms  with 
optional  choices  for  various  cases  has  been  presented  [14,15; 
and  will  be  published.   We  present  here  a  few  examples  of 
typical  formulations. 


A.   Axially  Symmetric  Resistive  Formulation  (original  Grad- 
Hogan) 

A*4;  =  -  r^dp/dij^  -  fdf/dii>    ,    f  =  rB 


||-   +  U-V4;  =  r]^*i)    -    c 


/^   ||N+  (f<u-VV/r^>)  •  -  (nKf)'   ,    •  =  3/3V 


p   +  Up'  +  YPU'  -    0      (Y=l  isothermal,  y>1  adiabatic) 


<A*iJj/r^>  -  (Ki|j')'   ,   K  =    <|VV|^/r^> 


Note:   The  averaged  9f/8t  equation  is  the  condition  for 

solvability  of  the  S-component  of  u  (which  component 
appears  only  in  the  equation  for  9f/3t) . 
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A, :   Formulas  with  V  as  Independent  Variable 

<l/r^>    =   A    ,    <r^><l/r^>    -    1    =    t^    ,    fA  =    x'     /    U   =    <u-VV> 
i>^   +    Uijj'    =     (n/A)  [{Ki>')  '    -    X^p'/H>'  ]    -    c 

Xt  +  ux'  -  nK(x'/A)'  -  nT^x'p'/A(ij;')^  -  c^ 

Two  interesting  identities  are: 

(19)  (x^  +  C^)/X'  -  i^^    +    c)/^'    -  nKf /fA  -  nCKij;')  '/Aij;' 

(20)  <B^>U  -  -  np'  (K<r^>  +  f^T^/(ij;')^)  -  Kip'  (ip^   +  c)  -  f  ( x^  +  c^) 

<B  >  =  <B   +  B^>   ,   <B  >  =  K<|;'     ,   <B^>  =  Af 

2 

With  "^      and  x^  dropped  (19)  represents  the  K  PS  profile 

restriction  and  with  the  time  derivatives  dropped  (20)  gives 
the  PS  transport  formula.   To  evaluate  the  global  GH 
"corrections"  requires  complete  solution  of  the  problem. 
This  formulation  is  incomplete  since  the  quantity  U  has  been 
neither  determined  nor  eliminated. 


A^ :  For  large  aspect,  small  poloidal  field,  finite  poloidal 
3  and  rotational  transform  (original  GH)  the  fast  time  scale 
evolves  according  to 


^^   =  (n/A)  (K^')  '  -  c 

and  all  terms  in  the  formula  for  U  (including  ^^   and  x^) 
contribute  equally  to  the  slow  pressure  diffusion  which 
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accompanies  the  fast  ^   diffusion.   After  ijj  equilibrates,  the 
i|;   and  x^    terms  drop  out  of  the  expression  for  U,  which  is 
essentially  Pf  irsch-Schliiter ,  and  gives  a  slow  equilibration 
for  p. 


A^:   For  finite  aspect  and  cross-section  and  finite  poloidal 
with  small  poloidal  field,  both  p  and  B  diffuse  at  the  same 
fast  rate, 


^^   =    (n/A)  (Kijj')  '  -  c 


2     2 

u  =  -  np'x  /Aijj' 


where  the  modified  classical  U  makes  p  diffuse  as  quickly  as 
ip.   The  two  diffusion  processes  seem  to  be  uncoupled,  but 
they  are  actually  coupled  by  the  pressure  balance-   These 
formulas  are  misleadingly  simple  because  the  coefficients 
K,  A,  T^  vary  on  the  same  time  scale  as  p  and  ijj  (unless  the 
cross-section  is  approximately  circular) . 


B:    Adiabatic  Axially  Symmetric  Formulation 
A*^  =  -  r^p  -  ff      ('  =  3/3i|^) 

p  =  u(t|;){4j')^   ,   f  =  v(4^)i|j'/A      (for  y-2) 

2  2 

(Ki|j')'  =  -  (y  +  vv/A)i|;'   -  2\iA^"    -    V  (tJ;'/A)' 

The  adiabatically  invariant  profiles  y(ij;)  and  v(ij;)  are 
obtained  by  eliminating  U  from  the  equations  for  p^,  4;^,  and  ^'^. 
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The  last  equation  (average  pressure  balance,  with  p  and  f 
removed  in  favor  of  M  and  v)  is  an  ordinary  differential 
equation  for  ^(V)    given  the  adiabatic  profile  y(i|)),  v{^)    and 
the  geometric  inductance  coefficients  K{V)  and  A(V).   The 
solution  procedure  involves  iteration  on  K{V)  and  A{V).   The 
2-D  elliptic  operator  is  inverted  with  known  right  side  after 
evaluating  a  profile  iJj  (V)  from  the  ODE  to  determine  a  revised 
set  of  contours  which  redetermine  K  and  A. 

Together  with  a  generalized  definition  of  the  concept  of 
adiabaticity  (allowing  islands  to  form  and  change  size,  subject 
to  appropriate  jump  conditions  across  a  Separatrix) ,  this 
algorithm  has  been  used  to  calculate  adiabatic  changes  of  a 
Tokaraak  and  adiabatic  transition  between  a  Belt  Pinch  and 
Doublet  [11].   With  added  heat  source  (-^  =  given),  it  has 

o  t 

also  been  used  to  describe  a  Flux  Conserving  Tokamak  [23] . 
The  pressure  balance  involves  all  relevant  inductance 
coefficients,  one  in  2-D,  two  in  axial  symmetry,  three  (full 
2^2  matrix)  in  helical  symmetry  [9] . 


C:    General  Tokamak,  Resistive  and  Heat  Conducting,  with 
Toroidal  Flux,  x  as  Independent  Variable 

To  the  i)      and  x^  equations  in  A,  we  add  the  energy  equation 

with  T  =  T{^) , 

p^  +  Up'  +  YpU'  =  (Y-1)(A^T')'  +  (Y-l)<nJ^> 

X^    =    A(T)<|Vv|2/aj2x2> 

<nJ^>    =    n[K(f')^    +    F^/A   +    T^(p')^/A(i(;')^]     ,    F   =     {Ki,')' 
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With  X  as  independent  variable, 


6fY  t)  •  ^1   =  ^    ^  = 


*t  +  U4)'  =  If  +  4'(Xt  +  ux') 


Let  x'  s  a/  y  =  P/a '  ,  \ 
As  an  abbreviation,  set 


i)    ,    E,    zz    p/a. 


1^. 


Ui|j' 


X 


up' 


YPU' 


Then, 


(C.I) 


dy  _ 
dt 

b^-^ 

dv 
dt 

(X    -    VY) 

dS    _ 
dt    - 

-     (CY)* 

YUY  -  yY 


The  pressure  balance  takes  the  form 


(C.2) 


av (Kav)   =  -  (ya  '  ) 


a (a/A) 


There  are  four  equations  for  (y,v,^,a)  ,  (y,v,C)  diffuse, 
while  a  satisfies  an  ODE.   The  derivatives  a  and  a  can  be 
explicitly  eliminated  in  favor  of  (y,v)  and  (y,v).   In  this 
way  the  second  derivatives  of  (y,v,5)  in  (C.l)  are  explicit, 
The  eigenvalues  of  the  3x3  matrix  of  second  derivatives 
determine  the  numerical  stability  and  accuracy  and  are  the 


•19- 


actual  diffusion  coefficients  for  small  wavelength  disturb- 
ances.  These  eigenvalues  (for  y  =  2)  are 

A^  =  nKa^/A 

where 

X^    =    (2na^u/AA) [K(l+T^)  +  x^/Av^] 

^2  =  (A^a/5) (1  -  y/A) 

6  =  ^   ,   A  =  Kv^  +  2y  +  1/A 

A   is  the  (geometry  corrected)  skin  diffusion  coefficient. 
There  are  many  interesting  consequences  in  the  limits  of  low 
and  high  3/  aspect,  etc.   Note  that  A   and  A_  are  respective- 
ly larger  and  smaller  than  the  larger  and  smaller  of  (A,,A-). 
Also 

A. 

"o 
A   is  the  plasma  density  diffusion  coefficient  if  the  heat 


^  =   6[1  +  T^(l+l/AKv^)]   ,   6  =  p/<|  B^> 


equation  is  replaced  by  a  prescribed  heat  source  dy/dt  =  Q, 
or  is  adiabatic.   In  general,  plasma  diffusion  and  heat  flow 
are  coupled  in  A_^.   Nonlinearities  couple  all  effects,  even 
when  the  eigenvalues  do  not  exhibit  this. 

This  code  has  been  programmed  for  a  simple  topology  [at 
ORNL]  and  is  being  programmed  for  a  Doublet  configuration  [at 
NYU]  . 
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D:    Resistive  Tokamak  (no  Heat  Flow)  with  Mass  as  Independent 
Variable 


The  notation  for  (j)(M,t)  ,  M  =  /pdx  is 


3(1)     _   d(j)      3(J) 
■^L,   ~   dt   '   9M 


H  ^  "*•  =  if 

Take  ^   and  x  as  diffusion  dependent  variables  and  p  as 
auxiliary  variable: 

II  =   (n/A)  [p(KpiI>)  *  -  T^p/h    -   c 


dt 


nKp(px/A)   -  HT  XP/Aip   -  c-|_ 


i\i{Kp^)       =  -  p/p  -  (x/A)  (PX/A) 

I,  p  =  y  (M)  p   where  y  (M)  is  ei 

satisfies  dy/dt  =  Q(M,t)  =  given.   The  third  equation  (ODE) 

can  be  used  to  eliminate  p  and  p  in  favor  of  ijj  and  x  to 

exhibit  these  second  derivatives  explicitly  in  the  diffusion 

equations  for  i|;  and  x*   The  eigenvalues  (diffusion  coefficients) 

are  A   and  \,    of  Sec.  C. 
o      1 

There  are  many  satisfactory  choices  of  variables,  but 
they  are  subject  to  stringent  restrictions.   In  order  of 
differentiation  (cf .  pressure  balance) ,  p  is  comparable  to 
ip'  and  x'«   Therefore  ii   and  x  can  be  taken  as  dependent 
variables  only  if  there  is  no  energy  equation  (case  D)  or  by 
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adjoining  /(p/P  ) dM  as  a  third  variable.   Or  (as  in  case  C) , 
moving  up  one  derivative,  (^^',x''P)  can  be  taken  -  but  only 
after  proper  normaliza-cion  (to  eliminate  U)  ,  viz.  \p'/p    i 
X'/P  I    p/P  •   In  normalizing  to  eliminate  U,  either  i|; '  or 


poloidal  field  configuration,  i|j  is  a  poor  independent  variable; 
in  a  simple  mirror,  x  -  0  cannot  be  used;  near  the  edge  of 

the  plasma  normalization  with  p,  e.g.  ij^'/p  is  unbounded, 

2 
while  p/p   has  no  evident  limit.   In  a  reversed  field  simple 

mirror  one  has  a  choice  between  ij^'/P  which  is  unbounded  at 

the  edge  and  p/i|^'  which  is  unbounded  at  the  separatrix.   The 

question  of  whether  a  boundary  condition  is  imposed,  or  left 

free,  or  whether  a  unique  natural  boundary  condition  must  be 

imposed  requires  individual  verification  for  each  variable 

at  the  plasma  edge,  center,  axis  (if  within  the  plasma),  and 

separatrix.   Many  (though  not  all)  of  these  permutations  have 

been  evaluated  and  will  be  reported  on  subsequently. 
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(a)  Initial  State 


ully  Developed 


Transition  From  Belt  Pinch  to  Doublet. 
Curves  of  Jvs.V  Show  Moving  Singularity, 


Figure  1 
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